When Does Division of Labor Lead to Increased System Output? 
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This paper develops a set of simplified dynamical models with which to explore the conditions 
under which division of labor leads to optimized system output, as measured by the rate of pro- 
duction of a given product. We consider two models: In the first model, we consider the flow of 
some resource into a compartment, and the conversion of this resource into some product. In the 
second model, we consider the resource-limited growth of autoreplicating systems. In this case, 
we divide the replication and metabolic tasks among different agents. The general features that 
emerge from our models is that division of labor is favored when the resource to agent ratio is at 
intermediate values, and when the time cost associated with transporting intermediate products is 
small compared to characteristic process times. We discuss the results of this paper in the context 
of simulations with digital life. We also argue that division of labor in the context of our replication 
model suggests an evolutionary basis for the emergence of the stem-cell-based tissue architecture in 
complex organisms. 

Keywords: Differentiation, division of labor, replication, metabolism, stem cells, tissue architecture, agent- 
based models 



I. INTRODUCTION 

Division of labor is a ubiquitous phenomenon in biol- 
ogy. In sufficiently complex multicellular organisms, var- 
ious tasks necessary for organismal survival (metabolism, 
nutrient transport, motion, reproduction, information 
processing,etc.) are performed by distinct parts of the 
organism 0, 0- Division of labor is even possible in 
clonal populations of free-living single-celled organisms 
At longer length scales, it is apparent that divi- 
sion of labor is a strong characteristic of community be- 
havior in various animals 0, 0. Human-built modern 
economies exhibit considerable division of labor (indeed, 
much research into this phenomenon has been done by 

economists) HHHUEflEJmm. 

Selective pressure for the division of labor in a popula- 
tion of agents (cells, organisms, humans) arises because 
specialization allows a given agent to optimize its per- 
formance of a relatively limited set of tasks. Total sys- 
tem production of a population of differentiated agents 
can therefore be significantly greater than a comparable 
population of undifferentiated agents. 

The question that arises then, is why is division of la- 
bor not always observed? For example, while complex 
multicellular organisms are certainly ubiquitous, approx- 
imately 80% of the biomass of the planet is in the form 
of bacteria. While capable of exhibiting cooperative be- 
havior, bacteria are, for the most part, free-living single- 
celled organisms. Clearly then, there are regimes where 
differentiation is not desirable. 

As a general rule, the more complex the organism, the 
greater the selective pressure for differentiation of system 
tasks. This rule is admittedly somewhat circular, since 
a more complex organism will by definition exhibit more 
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specialization of the component agents. So, to be more 
precise, the greater the number of agents comprising a 
system, the greater the selective pressure for differentia- 
tion (even this formulation has some ambiguity, because 
we can arbitrarily define any group of agents to com- 
prise a system, no matter how weak the inter-agent in- 
teractions. Nevertheless, despite this ambiguity, we will 
proceed with this initial "working" rule). 

The origin of this rule comes from the observation that 
there is a cost associated with differentiation, namely a 
time (and energy, though this will be ignored in this pa- 
per) cost associated with transporting intermediate prod- 
ucts from one part of the system to another. As sys- 
tem size grows, then presumably the density of agents 
grows (since the number of agents grows, and since we are 
grouping all the agents into one system, the inter-agent 
interactions are sufficiently strong, compared to some ref- 
erence interaction, to warrant this grouping. Note that 
increasing the agent density is a simple way to do this, 
though highly interconnected systems may interact fairly 
strongly over relatively long distances. The internet is a 
good example of this) . As the density of agents grows, the 
characteristic time associated with transporting interme- 
diates from one part of the system to another decreases, 
and so the cost of differentiation decreases (in fairness, 
the idea of transport costs placing a barrier to differentia- 
tion is not originally the author's. In the context of firms, 
this idea has been presented in the economics literature 

In this paper, we develop two sets of models that cap- 
ture the competition between the benefits of differentia- 
tion and the time cost associated with differentiation. In 
the first model, we consider the flow of some resource into 
a compartment, and its conversion into some final prod- 
uct. In the undifferentiated case, we assume that there 
is a single agent capable of converting the resource into 
final product. In the differentiated case, we assume that 
the conversion of resource is accomplished in a two-step 
process, each of which is carried out by distinct agents 
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specialized for the separate tasks. 

In the second model, we consider the flow of resource 
into a region containing replicating agents. We assume 
that the agents increase their volume so as to maintain a 
constant, pre-specificd population density. In the undif- 
ferentiated case, we assume that a given agent can absorb 
the resource, and process it to produce a new agent. In 
the differentiated case, we assume a division of labor be- 
tween replication and metabolism steps. That is, we as- 
sume that a fraction of the agents are specialized so that 
they cannot replicate, but can only process the resource 
into an intermediate form. This metabolized resource 
is then processed by the replicators, which produce new 
agents. These daughter agents then undergo a differen- 
tiation step, where they can either become replicators 
themselves, or metabolizers. 

In both the compartment and the replicator models, 
the general features that emerge is that differentiation is 
favored when population density is at intermediate levels 
with respect to resource numbers (when the population 
density is low, then the undifferentiated pathways are fa- 
vored, while when resources are highly limited, then the 
difference between the undifferentiated and differentiated 
pathways disappears). In the context of the replicator- 
metabolism model, we argue that this phenomenon sug- 
gests an evolutionary basis for the stem-cell-based tissue 
architecture in complex vertebrate organisms. 

This paper is organized as follows: In Section II, we 
develop and discuss the compartment model, involving 
the conversion of some resource into a final product. 
In Section III, we develop and discuss the replication- 
metabolism model. In Section IV, we conclude with a 
summary of our results and plans for future research. 

II. COMPARTMENT MODEL 
A. Definition of the model 

The compartment model is defined as follows: Some 
resource, denoted R, flows into a compartment of vol- 
ume V at a rate f R . In the undifferentiated case, a 
single agent, denoted E, processes the resource to pro- 
duce a final product, denoted P (the term E comes from 
chemistry, since the chemical analogue of an agent is an 
Enzyme catalyst). 

In the differentiated case, the processing of R is accom- 
plished by two separate agents, E\ and Ei. The agent E\ 
first converts the resource into an intermediate product 
R* , and then the agent Ei converts R* into P. 

It should be apparent that separating the tasks associ- 
ated with converting R to P among two different agents 
can only increase the total production rate of P if E\ 
and Ei can each perform their individual tasks better 
than E. Therefore, an implicit assumption here is that, 
when an agent specializes, its "native" ability to perform 
a given task can be made better than when an agent is 
unspecialized. 



For a simple reason why this is true, let us imagine 
that E, E\, Ei are enzymes, i.e. protein catalysts, whose 
function is pre-determined by some amino acid sequence 
of length L. If the alphabet is of size S, then there are 
S L distinct sequences that can generate E, E\, E 2 . As- 
suming that E\ and Ei are optimized for their particular 
functions, we note that, in the absence of any additional 
information, the probability that E\ and Ei are the same 
is 1/S L — > as L — > oo. Indeed, the average Hamming 
distance (number of sites where two sequences differ) be- 
tween any two sequences in the sequence space is given 
by L(l — l/S*) oo as L ^ oo. 

Therefore, it is highly likely that E is neither optimized 
for any of the tasks associated with converting R to P, 
but performs each task with some intermediate efficiency. 

B. Undifferentiated model 

In order to describe the processes governing the con- 
version of R to P, we will adopt the language and no- 
tation of chemical reaction kinetics. This formalism is 
very convenient, and is easily translatable into a system 
of ordinary differential equations. 

For the undifferentiated model, we have, 

E + R — > E — R second-order rate constant ki 

E — R — > E + P first-order rate constant k 2 

R —> Decay products (first-order rate constant kj$) 

The first reaction refers to agent E grabbing the resource 
R (in chemistry, this is referred to as the binding step). 
At this point, the agent is denoted E~R, to indicate that 
it is bound to a resource particle. In the second reaction, 
the agent processes the resource to form the product P, 
which it then releases. 

The last reaction indicates that the resource R has a 
finite lifetime inside the compartment, and decays with 
some first-order rate constant k D . This assumption en- 
sures that the compartment cannot be filled with resource 
without limit. The finite lifetime can be due to back- 
diffusion of resource outside the compartment, or simply 
that the resource does not last forever (some analogies 
include waiting time of customers at a restaurant before 
leaving without being served, the characteristic time for 
a food product to spoil, or the diffusion of solute out of 
a cell). 

If rift, ^e, tier denote the number of particles of re- 
source R, unbound agents E, and agent-resource com- 
plexes E — R, respectively, then we have, 

dn R ki 

—jj- = JR - -^n E n R - k D n R 

dn E ki 

—7— = --rrn E n R + k 2 n ER 
dt V 

dn ER h 

— 7— = T7 n EnR ~ km E R (2) 
at V 

If we define n — he + tieu, then note that dn/dt = 
dnE/dt + dnE R /dt = 0, which implies that n is constant. 
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After some manipulation, we obtain that the steady-state 
solution of this model is given by, 



nR, ss = 
where uer,ss satisfies 



k2riER,.s 



(ki/V)(n - n ER . ss ) 



(3) 



)n E R,ss + -i-n = (4) 

k 2 



so that, 



l r/ . Ir , k D f R k D 

nER > ss = ^VWJvj^f 71 + T 2 + WJvJ - 

(5) 

We take the "-" root because it guarantees that her, S s < 
n for all positive values of Jr. 

For small n, a Taylor expansion of the quadratic to 
first order gives, 



nER,ss — 



f R /k 2 



f R /k 2 + k D /(ki/V) 



n (small n) (6) 



while for large n, Taylor expansion to first order with 
respect to the remaining terms gives, 



f R 

n E R,ss = -r- (large n) 
k 2 



Note that the intermediate product, R*, is also capa- 
ble of decaying. As we will see shortly, it is the finite 
lifetime of the intermediate products that causes the un- 
differentiated pathway to outperform the differentiated 
pathway at low agent numbers, and allows for a transi- 
tion at higher agent numbers, whereby the differentiated 
pathway overtakes the undifferentiated pathway. 

In the context of our model, the direct interpretation 
of the decay term for R* is that the intermediate prod- 
uct has a finite lifetime, due either to diffusion out of 
the compartment or due to decay into other compounds. 
More generally, though, this term may refer to an aging 

cost, and therefore this model may be useful in under- 

4Z5^anding aspects of networked systems, whose function 
fedoes not necessarily depend on material transfers, but 
on information transfers. 

Information is transmitted between various parts of a 
system in order to effect system behavior, in response 
to the state of the system at the time of information 
transfer. Therefore, there is a time limit during which the 
information is relevant (because of the dynamic nature 
of the system and environment), which may be roughly 
modelled by assuming that the information is "lost" via 
a first-order decay. 

Defining particle and agent numbers analogously to the 
^ undifferentiated case, we obtain the system of equations, 



As a rough estimate of where the transition from the 
small n to the large n behavior occurs, we can equate the 
two expressions and solve for n. The result is, 



^trans.l , 



[r 
k 2 



CD 



(ki/V) 



(8) 



C. Differentiated model 

The conversion of resource R into P via the differen- 
tiated pathway occurs via the following sets of chemical 
reactions: 

Ei + R — > Ei — R second-order rate constant k[ 
Ei — R — > Ei + R* first-order rate constant k' 2 



E 2 
E 2 
R- 
R* 



R* 
R* 



E 2 
E 2 



R* second-order rate constant k' 3 
P first-order rate constant k\ 
Decay products first-order rate constant kn 
- Decay products first-order rate constant fcfp) 



dn R k[ 

—jj- = JR - -yn El n R - k D n R 

dn R* , k ' 3 
—— = k 2 n El R - -ryn E2 nR, - k D n R , 
at V 

dn El k[ , 
77 = ~~TF n Ei n R + k 2 n ElR 
at V 



dn El R _ k[ 

dt 
dn E2 



V 



n El riR - k 2 nE 1 R 



k' 

-yn E2 n R * + k' 4 riE 2 R* 
dn E2 R* k' 3 

— -j t — = -^n E2 n R , - k A n E2 R- 



(10) 



If wc define n\ — n El + nE 1 R and n 2 — ue 2 + riE 2 R*, 
then note that dni/dt = dn 2 /dt = 0, so that ni and n 2 
are also constant. Proceeding to solve for the steady- 
state of this model, we obtain, 
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ini 



;[(«2 



CD 



4' + /V) 



)-\(m 



k2"riE 1 R 
K 



+ 



(k'jvy 



>2 



(fci/V)' 



. /ft i 

rt9 



k' 2 n ElR 
K 



+ 



k* 



(11) 



r 



Now, when m and n 2 are small, it may be shown that 
to lowest non-vanishing order, the steady state popula- 
tion of E12 — R* is given by, 



while the small n expression for the rate of production of 
final product in the differentiated case is, 



nE 2 R*,ss 



{k'JV) k' 2 



f R /k' 2 



k* D K fn/k'2 + k D /(k[/V) 



a(l— a)n 2 (small n) 



k' 4 iiE 2 



k 



, (k' 3 /V) 



b* 



f R /k 2 + k D /(k[/V) 



a(l — a)n 2 



(12) 



where n = n\ + n 2 , and a = ni/n, and 1 — a = n 2 /n. 
For large values of n, we obtain that, 



n E2 R*, ss = ( lar S e n ) 



(13) 



As an estimate of where the transition between the 
small n and large n behavior occurs, we can equate the 
two expressions and solve for n. We obtain, 



1 



k 



Tltrans,2 



y 1,1 ^ 



a{l-a){k'jvyk' 2 (k[/V) 



) (14) 



D. Comparison of undifferentiated and 
differentiated models 

The small n expression for the rate of production of 
final product in the undifferentiated case is, 



k2TlER = k 2 



f R /k 2 



f R /k 2 + k D /(ki/V) 



(15) 



(16) 

Note then that for sufficiently small n, the undifferen- 
tiated production pathway produces final product more 
quickly than the differentiated pathway. However, be- 
cause the rate of production of final product for the un- 
differentiated pathway initially increases linearly with n, 
while the rate of production of final product for the dif- 
ferentiated pathway increases quadratically, it is possible 
that the differentiated pathway eventually overtakes the 
undifferentiated pathway. The critical n where this is es- 
timated to occur, denoted n equa i, may be estimated by 
equating the two expressions. The final result is, 

_ _J_ k* D fR/k'+kn/jkj/V) 

nequal ~ a(l - a) (k' 3 /V) f R /k 2 + k D /(h/V) U '> 
Now, for n equa i to be meaningful, it must occur in 
a regime where the rate expressions used to obtain it 
are valid. Therefore, we want n equal < n tranStl ,n tranSt2 . 
However, we can make an even stronger statement. If 
n eq uai does indeed refer to a point beyond which the dif- 
ferentiated pathway overtakes the undifferentiated path- 
way, then we should have n equa i < n t 

rans,2 < n trans,l- It 

is possible to show that, 



Tltrans,2 Tltrans,l 

T^equal Tltrans,2 



= Ja(l-a) 



(k' 3 /V) 



k D 



k* D f R /k' 2 + k D /{k'JV) y k 2 (h/V) 



;) 



(18) 



r 



and so, our inequality is equivalent to the condition that, 

k D J R /k2 + k D /(k 1 /V) 



kh <a(l-a)( fR 



(k's/V) 



k 2 {h/vy f R /k' 2 + k D /(k[/v) 

(19) 



f R 



which implies that 

n equa i<~-^ , 
Figures 4 and 5 (of the version submitted to The Jour- 



(20) 



nal of Theoretical Biology) show comparisons of the pro- 
duction rates of final product for the undifferentiated and 
differentiated pathways. In Figure 4, the parameters are 
chosen so that the differentiated pathway eventually over- 
takes the undifferentiated pathway, while in Figure 5, the 
parameters are chosen so that this is not the case. Note, 
however, that even in Figure 4, although the differen- 
tiated pathway overtakes the undifferentiated pathway, 
once n becomes very large, the undifferentiated pathway 
again overtakes the differentiated pathway. 
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This behavior can be explained as follows: When n is 
very small, the rate at which the intermediate product is 
"grabbed" by E 2 is small compared to the decay rate, so 
that much intermediate product is lost. In this regime, 
the undifferentiated pathway is optimal, for, although 
E may be less efficient than either E\ or E 2 at their 
respective tasks, the overall production rate of P is not 
reduced by the loss of intermediates. 

Now, as n increases, the rate of loss of intermediates 
decreases to an extent such that the increased efficiency 
associated with differentiation causes the differentiated 
pathway to overtake the undifferentiated pathway. How- 
ever, once n increases even further, then there is sufficient 
quantity of agents in both the undifferentiated and differ- 
entiated pathways to process all of the incoming resource, 
with minimal loss due to decay. At this point, because 
the production rate of P has become resource limited, 
the efficiency advantage of the differentiated pathway is 
considerably reduced, such that the slight cost associated 
with intermediate decay becomes sufficient to cause the 
undifferentiated pathway to overtake the differentiated 
pathway. However, this effect is a small one, since, once 
n is very large, both the undifferentiated and differenti- 
ated pathways perform similarly. 



E. When can a differentiated pathway outperform 
an undifferentiated pathway? 

The analysis of the previous section deserves further 
scrutiny, in order to better understand the circumstances 
under which differentiation can lead to improved system 
performance. 

At low agent numbers, the decay of the product inter- 
mediates leads to a quadratic increase in system output, 
and so the undifferentiated pathway outperforms the dif- 
ferentiated pathway. At some point, however, the num- 
ber of agents is sufficiently large that the decay of both 
resource and intermediates is minimal, so that it is pos- 
sible for the differentiated pathway to overtake the un- 
differentiated pathway. However, this is only possible if 
the differentiated pathway is more efficient than the un- 
differentiated pathway. 

To quantify this notion, assume that n is at some in- 
termediate value, such that ku and k* D may be effectively 
taken to be 0. In this regime, it is possible to show that, 



k' 4 nE 2 R* , S s = min{fc 2 cm, k' 4 (l — a)n} (21) 

Essentially, if k' 2 an > k' 4 (l — a)n, then the first set of 
agents are capable of producing intermediate at a rate 
greater than the second set of agents are capable of pro- 
cessing it, so that the second reaction step is rate limiting. 
If k' 2 an < k' 4 (l — a)n, then the first reaction step is rate 
limiting. 

Note then that if one of the reactions is rate limiting, 
we can adjust the agent fractions to increase the rate 
of the rate limiting reaction, and thereby increase the 



overall production rate of P. Therefore, the maximal 
production rate of P is achieved when k' 2 an = k' 4 (l — 
a)n =4> ot op timai — k' 4 /(k 2 + k' 4 ), so that the maximal 
production rate of P is given by, 



(k 4 n E ., 



2-R* ,ss jvaax 



k 2 k 4 



k' 2 



(22) 



For the undifferentiated case, the analogous expression 
is k 2 n, and so we expect that the differentiated pathway 
can only overtake the undifferentiated pathway when, 



K 2 K 4 

k 2 -\- k 4 
I 

k' 2 



> k 2 



k 4 



< 



1 

~k 2 



(23) 



Intuitively, this condition makes sense, since l/k 2 is 
the characteristic time it takes agent E to convert R 
to P, while l/k' 2 and l/k' 4 are the characteristic times 
for agents E\ and E 2 to perform their respective tasks. 
Therefore, differentiation can only overtake nondiffcrcn- 
tiation if the characteristic time for the completion of a 
set of tasks is shorter for the differentiated pathway than 
it is for the undifferentiated pathway. 

If 1/&2 + > it is in principle possible for 

the differentiated pathway to overtake the undifferenti- 
ated pathway if k[ is sufficiently greater than k\, and if 
k' 3 is sufficiently large compared to k* D . Basically, the 
differentiated agents are not more efficient at actually 
processing the resource, but they are more efficient at 
grabbing them, which can give the differentiated path- 
way an advantage. However, in contrast to the case 
where l/k' 2 + l/k' 4 < l/k 2 , this advantage is only tempo- 
rary, because once the agent number becomes sufficiently 
large, the characteristic time to grab resource becomes 
very small. 

As a final note, because the condition for differenti- 
ation to outperform nondifferentiation at larger agent 
numbers is l/k 2 + l/k' 4 < l/k 2 , while the agent number 
n eqU ai where differentiation overtakes nondifferentiation 
does not depend on k' 4 , it should be apparent that our 
criterion for n equa i could be inaccurate in actually pre- 
dicting the location of the cross-over. This is because 
n eqU ai is based on the small n region, where the rate of 
production of P for the differentiated pathway increases 
quadratically. This is the regime where the production 
rate of P is limiting by intermediate resource decay. 

However, as l/k' 2 + l/k' 4 increases, we expect n equa i to 
become a better predictor of the crossover point, since 
the decay of the intermediate resource becomes a com- 
paratively greater factor in dictating the performance of 
the differentiated pathway. 

In any event, though, the expression for n equa i and the 
condition for the existence of a cross-over are useful, for 
they indicate that the larger the value of Jr, the larger 
the value of k* D that is possible for a cross-over to still 
occur. In particular, it suggests that, as long as l/k' 2 + 
1/&4 < f/fc 2 , then by making fa sufficiently large for a 
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given k* D , we will eventually obtain that the differentiated 
pathway will outperform the undifferentiated pathway at 
sufficiently high agent numbers. This is indeed what is 
observed numerically. 



III. REPLICATION-METABOLISM MODEL 

In this section, we turn our attention to the replication- 
metabolism model, where a population of agents pro- 
cesses an external resource for the purposes of producing 
more agents. 



A. Definition of the model 

We consider a population of replicating agents, relying 
on the supply of some resource, denoted R. We assume 
that the resource is supplied to the population at a rate 
of f R per unit volume, and that, as the population grows, 
the volume expands in such a way as to maintain a con- 
stant population density p. 

In the undifferentiated model, a single agent, denoted 
E, processes the resource R and replicates. In the differ- 
entiated model, an agent E\ "metabolizes" the resource 
to some intermediate R* , and then another agent, de- 
noted E2, processes the intermediate and reproduces. 
However, the E 2 agents are responsible for supplying 
both metabolizers and replicators. Therefore, E2 pro- 
duces a "blank" agent, denoted E, which then specializes 
and becomes either E\ and E 2 . 



B. Undifferentiated model 

The reactions defining the undifferentiated model are, 

E + R — > E — R second-order rate constant k\ 

E — R — > E + E first-order rate constant k 2 

R — ► Decay products (first-order rate constant Hjify) 

In terms of population numbers, the dynamical equa- 
tions for tie, n ER and n R are, 

dn E fci 
at V 
dn E R h 

77 = -T7 n E n R ~ k 2 n ER 

at V 

= IrY - yn E n R - k D n R (25) 



Therefore, defining n = tie + n ER, we have, 
dn 

— = K 2 n E R 
at 



(26) 



Since the population density is p, this implies that, 
dV ldn 



dt p dt 



= k 2 Vx E R 



(27) 



where x E = n E /n, x E r = tier/ti. 

Now, the concentration c R of the resource R is given 
by the relation n R — crV, which implies that 

dc R 1 , dn R dV 
dt ~ V { dt CR dt> 

= !r - {kipx E + k 2 x E R + k D )c R (28) 

Putting everything together, we obtain, finally, the sys- 
tem of equations, 



1 dn 
n dt 
dxE 
~~dT ' 
dx E R 

dt 
dcR 
dt 



= k 2 x E R 

= ~kic R x E + 2k 2 x E R - k 2 x e rx e 

= k x c R x E - k 2 x E R - k 2 x\ R 
= Jr - c R {k x px E + k 2 x E R + k D ) (29) 



We can determine the steady-state behavior of the 
model by setting the left-hand-side of the above system 
of equations to 0. When p — > 0, the steady-state solution 
is characterized by, 



cr,. 



h 



(p-0) 



k D + k 2 x E R 
where xer,ss is the solution to the cubic, 



(30) 



kp ^ 2 , fkifn , k D ^ kif R 
k 2 " 



=a; 3 + ( 1 + ^)^ + (™ + ^l )a; _™ (/ ^ () , m] 



Note that when f R = 0, we obtain xer, ss — 0. Therefore, 
differentiating the cubic with respect to x gives, 



r dx E R, ss , ki 
-)f R =o 



and so we have, 

XER,s! 



dff 



ki 
k 2 k_ 



k 2 k D 



2«D 



f R (small f R , p^0) 



(32) 



(33) 



When f R is large, we get xer, ss — ► 1, so that 

x E R, ss = 1 (Jr. 00, p 0) (34) 

Equating the small f R and large f R expressions, we 
obtain that the transition from small f R to large f R be- 
havior is approximated by, 



fR,trans,l(P = °) = k 2 



D 



ki 



(35) 



Now, when p is large, then the steady-state expression 
for c R is approximated by, 

= ~JT = fR _ c R kipx E k x c R x E = — (36) 
dt p 
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and so, 



= - k 2 XER,ss - k2X 2 ER s 



(37) 



from which it follows that, 



2 y fc 2/ o 



(38) 



Since p is large, we will approximate this expression fur- 
ther, by taking the first-order expansion in f R / (k 2 p), giv- 
ing, 



XER.: 



k 2 p 



(39) 



We can estimate the cross-over from small p to large p 
behavior by equating the two expressions. We have two 
estimates, one for small f R , and one for large f R . We 
obtain, 



f- 'trans, 1 ~ (/« < 

Plans,! = ^ UR > *2&) 



(40) 



C. Differentiated model 

The reactions defining the differentiated model are, 

Ei + R — > Ei — R second-order rate constant k[ 
Ei — R — > Ei + R* first-order rate constant k 2 
E 2 + R* — > E 2 — R* second-order rate constant k' 3 
E 2 — R* — > E 2 + E first-order rate constant k' 4 
E — > Ei first-order rate constant k' 5 
E — > E 2 first-order rate constant k' 6 
R — > Decay products (first-order rate constant ko) 
R* — > Decay products (first-order rate constant k$$) 



Following a procedure similar to the one carried out 
for the undifferentiated model, we obtain the system of 



equations, 

1 dn 

n dt 
dx El 



dt 
dx El R 

dt 
dx E2 

dt 

dXE 2 R 



k' 4 XE 2 R* 

z -k[cRX El + k' 2 x El R + k' 5 XE - k' 4 XE 2 R*x El 
= k' 1 c R x El - k' 2 x El R - k 4 x E2 R*x El R 



= -k' z C R <>XE 2 + k' 4 X E2 R* + k' 6 XE - k' 4 XE 2 R*XE 2 

= ko > CR*XE 2 — k 4 XE 2 R* — k 4 x E2 R* 



dt 

—j^ = k' 4 XE 2 R* ~ (k' 5 + k' 6 )x E - k 4 x E2 R'X E 
—^j- = Jr. - (k'ipx El + k D )c R - k 4 XE 2 R*c R , 



dc R * 
dt 



p(k' 2 x El R - k' 3 c R *x E2 ) - (k* D + k 4 x E2 R*)cR. 



(42) 



Now, defining xe x — xe x + xe ± r, and xe 2 — xe 2 + 
xe 2 r*, we obtain, 



dx El 

dt 
dxE 2 

dt 



k' 5 XE - k' 4 XE 2 R*XE 1 
k' 6 XE - k' 4 XE 2 R*XE 2 



Therefore, at steady-state we have, 



XE u s 



4 



(43) 



(44) 



XE 2 ,ss n-6 

and so, using the relation x El + Xe 2 + Xe = 1 we obtain, 



x E lt 



xe 2 



7(1 - x Ey 



ka 



k' 5 + fcg 



7(1 - X E ,ss) 



(45) 



If we let k' 5 , fcg — > 00 such that k' 5 /k' 6 remains constant, 
then it should be clear that xe, ss —* 0. Intuitively, E dif- 
ferentiates to either E 4 or E 2 as soon as it is produced, 
so it does not build up in the system. The ratio between 
k' 5 and k' Q then dictates the fraction of Ei and E 2 in the 
system (allowing k' 5 , k' 6 — ► oo essentially amounts to as- 
suming that the differentiation time is zero. This is of 
course not true, and future research will need to incor- 
porate positive differentiation times). 

Defining a = k' 5 /(k' 5 + k' 6 ), we then have xe 1 , ss = a, 
&ndxE 2 ,ss = 1 — a. Therefore, to characterize the system 
at steady-state, we need to solve four equations, giving 
the steady-state conditions for x Ei r, xe 2 r* , cr, and cr- , 
respectively. The equations are, 

= k[c R (a - x Ei r) - k' 2 x El R - k' 4 x E2 R*x El R 
= k' 3 c R ,(l -a- x E2 r>) - k! 4 x E2 R- - k' 4 x% 2R , 
= /r - c R (k D + k[p(a - x Ei r)) - k' 4 x E2 R*CR 
= p(k' 2 x El R - k' 3 c R »(l - a- Xe 2 r*)) 
-k* D c R * - k' 4 c R ,x E2 R* (46) 
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As with the undifferentiated case, we study the behavior 
of this system of equations in both the small and large p 
limits. 

When p = 0, we have cr. iSS = => Xe 2 r*, ss = 
cr.ss = IrI^d => x Ei r. ss = {k' 1 f R a/k D )/(k' 2 + 

Kf R /k D ). 

Differentiating the steady-state equations with respect 
to p, and evaluating at p = 0, gives, 

(-^-)„=o - -gr (XE lR ) P =0 

,dxE 2 R',ss ■> k 3 d,CR* ss 

{ ^p— ) ^ = v} l - a){ ^p- )p=0 (47) 

and so, for small p, we have, 



XE 2 R*,ss — , , 



FT (small p) (48) 

/c 4 fc D k' 2 + ^ 

Now for large p, our steady-state equations may be 
reduced to, 

= k[c R (a - x Ei r) - k' 2 x El R - k' 4 x E2 R*x El R 
= k' 3 c R * (1 - a - xe 2 r*) - k' i XE 2 R* - k' 4 x E2R » 
= f R - k[cRp(a - x Ei r) 

= k 2 x El R-k 3 CR*(l-a-x E2 R*) (49) 

The third equation gives k[cR(a — Xe 1 r) — Ir/p, 
which may be substituted into the first equation to give, 



= — - x ElR {k' 2 + k' 4 x E2 R* 
P 



(50) 



Solving for Xe 1 r in terms of xe 2 r*, and plugging the re- 
sulting expression into the fourth steady-state equation 
gives, after some manipulation, that Xe 2 r*, S s is the solu- 
tion of the cubic, 



y 

n 3 /i 2 \ 2 2 

= x E 2 R\ss + (1+ -rr) x E 2 R%ss + -jjXE 2 R*,ss ^ 



^2 



(51) 

Now, when f R = 0, we obtain Xe 2 r*.ss — 0. From this 
it is possible to show that, 



and so, 



,dXE 2 R".ss^ 1 

[ df R )f "=° ~ ~pg 



xe 2 r*,ss = -n- ( lar S e P) 
k iP 



(52) 



(53) 



As with the undifferentiated case, the transition from 
small p to large p behavior may be estimated by equating 
the two expressions and solving for p. The result is, 



- 1 k D k h , 1 , fR n fr^ 



D. Comparison of undifferentiated and 
differentiated models 



As a function of we wish to determine if, as p 
increases, the average growth rate of the differentiated 
population overtakes the growth rate of the undifferenti- 
ated population. If this does indeed happen, then there 
exists a p, denoted p e quai, at which the two growth rates 
are equal. 

We first consider the regime fx < k^ (ku / k\ ) . This is 
the small Jr regime of the undifferentiated population. 
In this regime, the transition from low p to large p be- 
havior occurs at ptrans.i = ko/ki. For the differentiated 

pathway, we have p tr ans,2 = \J a (i- a) V^^Jr + k ~o)- 
Now, we have four possibilities: (1) p eqU ai < 

Ptrans,lj ptrans,2- (2) ptrans,2 < Pequal < Ptr 

(3) ptrans.i < Pequal < Ptrans.2- (4) Pequal > 
Ptrans.i? Ptrans,2- 

We can immediately eliminate Cases (2), (3), and 

(4) as possibilities. For Case (4), we get an undiffer- 
entiated rate of fe/flAfep) = f R /p, and a differenti- 
ated rate of k A fR/(k' 4 p) — Jr/p, and so the two rates 
are equal. For Case (2), we get an undifferentiated 
rate of kif R /k E , and a differentiated rate of Jr/p, so 
equating gives p equa i = k D /k x = p tra ns,i- Therefore, 
Case (2) is essentially a limiting case of Case (4), and 
can also eliminated. For Case (3), we get an undif- 
ferentiated rate of Jr/p, and a differentiated rate of 
(k 2 k 3 /k* D )(k' 1 f R /k D )/(k' 2 + k' 1 f R /k D )a(l-a)p, so equat- 
ing gives pequal = Ptrans.2- Therefore, Case (3) is essen- 
tially a limiting case of Case (1), and can also be elimi- 
nated. 

For Case (1), we have, 



Pequal 



1 hk* , 1 



Jr 



a(l — a 
Now, we can show that, 



fc 3 



(55) 



Pequal Ptrans,2 



k\\ 



1 



\ 1,1 U k„ ' 



Ptrans.2 Ptrans,! \j Ol(l - a) k' 3 k D k[ k' 2 k D ' 

(56) 

and so, in order for p eqU ai < Pt rans.li Ptrans,2i then we 
must have, 



k 2 k 2 r . * 

Jr < k D — — [a(l - a) 
k\ k 2 



k 3 kp 
k\kj-) 



(57) 



We now consider the case where Jr > korr- This is 
the large jr regime of the undifferentiated population. 
Following a similar procedure to the one carried out for 
the small Jr regime, we can show that the only possible 
crossover occurs in the small p regimes for both the un- 
differentiated and differentiated cases. In this regime, we 
obtain, 



Pequal 



1 



k D k* D A 



h 



f R a(l - a) k' 3 k[ k' 2 k D 



) 



(58) 
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We can show that, 



Pequal Ptrans^l 

Ptrans,2 Ptr ans.l 



= frack 2 fR\ 



a(l — a) 



k D k* D 1 



and so, in order for 
have, 



Pequal ^ Ptrans^li Ptrans,2i 



k 3 



k D /K + f R /k' 2 



V 

K \ "-2 

(59) 
we must 



(60) 



the availability of resource: When resources are plenti- 

fail, then increasing the resource favors a differentiated 

Jr replication strategy. However, when resource are lim- 
fco&i^ed, then decreasing the resource favors a differentiated 



In Figure 9, we show a high-//} plot where the differen- 
tiated growth rate overtakes the undifferentiated growth 
rate. In Figure 10, we show a high-/jj plot where the un- 
differentiated growth rate stays above the differentiated 
rate at all values of p (these figures are included in the 
version submitted to The Journal of Theoretical Biology) . 



E. When can a differentiated population 
outreplicate an undifferentiated population? 

We can subject our replication-metabolism model to a 
similar analysis to the one applied to the compartment 
model. First of all, as with the compartment model, 
we expect that the differentiated pathway can only over- 
take the undifferentiated pathway, and then maintain a 
higher replication rate if l/k 2 + l/k' 4 < l/fc2- Again, this 
condition simply states that the total characteristic time 
associated with converting resource into a new agent in 
the differentiated case is less than the total characteristic 
time in the undifferentiated case. The assumption is that 
decay costs are negligible, as well as time costs associated 
with grabbing resource and intermediates. 

An interesting behavior that occurs with the 
replication-metabolism model is the different dependence 
on Jr that the transition population density p eqU ai has 
in the Iow-Jr regime and the high-/;? regime. 

In the high-fR regime, p equa i has a weak dependence 
on Jr, though it does decrease as /r increases. This 
makes sense, for in the high- /a regime, the growth rate 
of the undifferentiated population is limited by the rate 
at which the complex E — R produces new agents. As Jr 
increases, the cost associated with the decay of the inter- 
mediate resource R* decreases, so that the differentiated 
pathway overtakes the undifferentiated pathway sooner. 

In the Iow-Jr regime, p equa i increases linearly with Jr, 
so that, as Jr increases in this regime, the differentiated 
pathway overtakes the undifferentiated pathway only at 
higher values of p (if it overtakes at all). The reason for 
this behavior is that at low Jr, the growth rate of the 
undifferentiated population is resource limited, so that 
increasing /r actually increases the growth rate. The ef- 
fect of this is to push to higher values of p the point at 
which the differentiated agents outreplicate the undiffer- 
entiated agents. 

What is interesting with these patterns of behavior 
is that they indicate opposite criteria for when a co- 
operative replicative strategy is favored, depending on 



replication strategy. 

In this vein, it is interesting to note that complex mul- 
ticellular life is only possible in relatively resource-rich 
environments. On the other hand, organisms such as 
the cellular slime mold (Dictyostelium discoideum) tran- 
sition from a single-celled to a multi-celled life cycle when 
starved. While we have already postulated one possible 
reason for this behavior in terms of minimizing overall re- 
productive costs Q , the behavior indicated in our model 
may provide another, complementary explanation for the 
selective advantage for this phenomenon. 



IV. CONCLUSIONS AND FUTURE RESEARCH 

This paper developed two models with which to com- 
pare the performance of undifferentiated and differenti- 
ated pathways. The first model considered the flow of 
resource into a compartment filled with a fixed number 
of agents, whose collective task is to convert the resource 
into some final product. The second model considered 
the replication rate of a collection of agents, driven by 
an externally supplied resource. 

By assuming that the resource, and even more impor- 
tantly, that reaction intermediates, have a finite lifetime, 
we were able to show that undifferentiated pathways are 
favored at low agent numbers and/or densities, while dif- 
ferentiated pathways are favored at higher agent num- 
bers and/or densities. An equivalent way of stating this 
is that differentiation is favored when resources are lim- 
ited, where resource limitation is measured by the ratio 
of available resource to agents. 

Some interesting results that emerged from our studies 
was that, although limited resources favor differentiation 
(as measured by the resource-agent ratio) , for a given set 
of system parameters, differentiation will be more likely 
to overtake nondifferentiation at higher population size 
and/or density if the amount of available resource is in- 
creased (although the actual cross-over location will in- 
crease as well). The central reason for this is that the 
relative decay costs associated with differentiation are 
decreased as resource is increased. 

In the context of the replication- metabolism model, we 
should note that when resources are plentiful, differenti- 
ation is favored at lower population densities as the re- 
source flow is increased, while when resources are limited, 
differentiation is favored at lower population densities as 
the resource flow is decreased. Regarding the former ob- 
servation, it should be noted that it has been shown that 
diversity of replicative strategies is favored at interme- 
diate levels of resources In digital life simulations, 
they showed that the number of distinct replicating com- 
puter programs was maximized at intermediate resource 
availability, a result consistent with what is observed eco- 
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logically. We claim that the results of this paper are 
consistent with these observations. 

Regarding the latter observation, we pointed out in the 
previous section that this behavior is possibly consistent 
with the behavior of organisms such as the cellular slime 
mold, which transition from a single-celled to a multi- 
celled life form when starved. 

We also posit that the results of the replication- 
metabolism model suggest a possible evolutionary ba- 
sis for the stem-cell-based tissue architecture in complex 
multicellular organisms. Essentially, as population den- 
sity increases, and therefore as the resource-to-agent ra- 
tio decreases, it becomes more efficient for some cells to 
exclusively focus on replicating and renewing the popu- 
lation, while other cells engage in specialized functions 
necessary for organismal survival. 

Of course, our replication-metabolism model is not 
quite the same as a stem-cell-based tissue architecture. 
First of all, the stem-cell and tissue cell population does 
not collectively grow. Rather, the stem cells periodically 
divide in order to replace dead tissue cells. Therefore, 
the stem-cell-based tissue architecture is a kind of hybrid 
between our compartment model and our replication- 
metabolism model. 



Secondly, our replication-metabolism model assumes 
that there is a single differentiation step, while in reality 
a differentiating tissue cell undergoes several divisions 
and differentiation steps before becoming a mature tissue 
cells. 

Finally, our replication-metabolism model assumed 
that differentiation was instantaneous. In reality, dif- 
ferentiation takes time, and this time cost will affect 
whether differentiation can overtake non-differentiation, 
and, if so, will likely delay the critical population density 
where this happens. 

Despite these shortcomings, we believe that the models 
developed here could be used as the basis for more sophis- 
ticated models that could produce, via an optimization 
criterion, the stem-cell-based tissue architecture observed 
in complex multicellular organisms. This is a subject we 
leave for future work. 
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